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Abstract 


In the present work, first, a procedure has been developed for the simulation of 
viscous diffusion in a boundary-layer flow. An operator splitting algorithm has 
been developed in which the viscous diffusion and convection are treated one after 
the other. The random walk technique is used for simulation of viscous diffusion. 
The boundary-layer velocity profile is calculated at one station and compared with 
theoretically obtained profile. 

The second part of the thesis is the simulation of the flow past a bluff body using 
the full-vortex cloud modelling by the surface vorticity boundary integral method. 
It develops a discrete vortex element model wherein the body is represented by a 
number of discrete vortex panels. The strengths of the vortices are related directly 
to the tangential velocity outside the boundary-layer. Finally the drag and lift 
coefficients are calculated by getting pressure distribution on the surface of the body 
and compared with the available data. We get good representation of the nature of 
the wake but the pressure calculation suffer from too much numerical noise. 

In the last part of this thesis an attempt is made to develop a procedure for cal- 
culating the flow past multiple two-dimensional bodies. The assembly technique of 
two mutually interacting bodies is found to be encouraging. We got the satisfactory 
vortex cloud solution for various Reynolds numbers. The main disadvantage in this 
method is that there is a lot of numerical noise in pressure calculation in unsteady 
flow, and the results are not yet satisfactory. 
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Chapter 1 


INTRODUCTION 


1.1 Literature survey 

The postulate by L.Prandtl in 1903 that fluid flows at large values of Reynolds 
number can be modelled as essentially inviscid with thin regions of viscous activities 
near the solid boundaries could not be used with profit in flow about bluff bodies, 
where the separation of these thin regions of viscous flow, termed as boundary 
layers from the solid surfaces disturbed the flow rather drastically. Consequently, 
the inviscid methods which depended upon the fact that regions of viscous flows were 
thin and close to the body and the pressure variations across them was negligible, 
could not be of much use in the calculation of flow about the bluff bodies. The 
discrete- vortex methods are a class of methods that retain the chain and simplicity 
of the inviscid method, and yet are able to account for the development of wake 
behind bluff bodies. The discrete-vortex methods have evolved slowly and have 
built upon the concepts introduced by various researchers. 

Theodore Von Karman had shown in 1911 that at relatively low values of Reynolds 
numbers the wake behind a bluff body consists of a series of vortices, which are 
shed alternately from the two separation points, one on either side of the body, and 
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convected downstream with the fluid. It was natural to explore the dynamics of such 
vortices and to exploit it for predicting the flows about bluff bodies. H elmh oltz has 
shown earlier in 1858 that in an inviscid flow, vortex lines are material lines, i.e. they 
are continually composed of the same fluid particles, and that flows with vorticity can 
be modelled with vortices of appropriate circulation and in fini tes im al cross-section. 
It was this realization which later led to the discretization of compact regions of 
vorticity into a number of lines or points vortices embedded within potential flows. 

Rosenhead (1931) is generally recognized to be the first person to use this 
approach to study the problem of Kelvin-Helmholtz instability of vortex sheets. 
He discretized a vortex sheet into a number of discrete vortices and allowed- these 
vortices to convect under mutual influence after an initial sinusoidal perturbations. 
After a number of time steps, the vortex roll up along the sheet developed, culminat- 
ing finally in a simple row of vortex bunches, very similar to half of a von Karman 
vortex sheet. 

After the success of Rosenhead, Abernathy and Kronauer (1962) extended this 
method to modelling the alternating vortex shedding in the wake of the circular 
cylinder by considering the stability of two parallel vortex sheets. Under the self- 
convection superposed on a free-stream velocity, the Karman sheet like pattern soon 
evolves. Abernathy and Kronauer concluded that the role of the bluff body is not 
very important in the formation of the characteristic vortex street except in its 
function as a source of the vortex sheet. 

There has been a large number of studies which have used various schemes to 
solve a large variety of problems. The excellent reviews by Clements and Maull(1975), 
Saffman and Baker (1979), Saffman (1981), Leonard (1980) and Sarpakaya (1989) 
summarized the activities of several groups. Much of the physical insights into the 
phenomenon of vortex shedding that has formed the basis of the vortex methods 
have been reviewed by Windhall (1975) and Mair and Maull (1971). 
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1.2 Vortex cloud approach 

The method of vortices consists of introducing in an essentially potential isolated 
vorticity elements (of various descriptions) which then convect under the superposed 
potential flow and the mutual influence of the vortices which can be predicted by 
Biot - Savart law or Green’s function method. The vorticity elements are tracked 
numerically using Lagrangian (or mixed Lagrangian - Eulerian) Schemes. 

1.3 Shedding of vortices 

The application of discrete vortex method to separated flow about bluff bodies 
requires that modelling of the separating vortex sheet, and its discretization. In 
general, two classes of methods can be identified those involving separation at fixed 
position in which the location of separation points is allowed to be determined by 
some condition. 

Clement (1973) and Kiya and Arie (1977, 1988) preferred the use of fixed nasent 
vortex position, where new vortices appear periodically in time at fixed spatial 
locations close to the experimentally known separation points. The rate of vorticity 
generation is obtained from the shear layer velocity U s at these points. 

§ = 0.5t? (1-1) 

This formulation is sought to be justified from the pioneering work of Fage and 
Johansen (1928) on structure of vortex sheets. Kiya and Arie (1977) carried out a 
parameteric evaluation of the variations of the location of these nasent - vortex loca- 
tion points and compared the calculated results with those obtained experimentally 
to fix a proper location. 

Sarpakaya (1989) refers to a work of Kuwahara (1973) wherein he determines 
the strengths of the nascent vortices appearing at fixed locations by using the Kutta 
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conditions. Kuwahara is reported to have found that the amplitude of the normal 
force coefficient calculated in this manner depends very strongly upon the assumed 
position of the two fixed points where the nascent vortices appear. In fact, a small 
change in the position of nascent vortices is reported to cause an eight - fold change 
in the normal force coefficient. 

In these methods, no interaction is allowed between the shed vortices and the 
appearance of the nasent vortices. Thus, the velocities at the outer edges of the 
shear layer are only indirectly related to the strength of the vortex and the fixed 
time interval. Evidently, the velocities in the inner and the outer edges of the shear 
layer, the strength and position of the nascent vortices, and the Kutta condition 
are interdependent and, therefore, the location and time interval cannot be chosen 
arbitrary. 

Sarpakaya (1975) introduced a procedure wherein he used variable nascent-vortex 
positions using Kutta condition. Parag Kumar (1993), using the lead provided by 
Sarpakaya based the strength of the nascent vortices on 

§ = 0.5 Ul„ (1.2) 

where U SiTn is the mean velocity of the shear layer calculated using the average of 
the convective velocities of the first four vortices in each shear layer. This appears 
to provide a satisfactory mechanism for feedback to the rate at which the vorticity 
is shed. 

1.4 Body Representation 

For calculating flows past bodies, the vortex methods, like all other inviscid methods, 
require that the flow tangency condition be satisfied. The no-slip condition is another 
matter and will be discussed later. 
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For many 2-D geometries conformal transformation offers a convenient method 
for satisfying the tangency condition. For a circular cylinder a simple transformation 
and the use of circle theorem yields one image per vortex at the inverse point. For a 
cylinder of arbitrary cross-section Sarpakaya (1975), Clement and Maul (1975) used 
Shwartz- Christoffel transformation, Shoff and Franks (1981) used truncated Laurent 
series, and Parag Kumar(1993) used a Zukowsky transformation for calculating flow 
past a flat plate at various angles of attack. 

The use of conformal transformation for 2-D geometries involve repeated trans- 
formation from physical plane to circular plane and to W-plane for velocities must 
be calculated in the W-plane but the convection of vortices must be calculated 
in the physical plane. Further Routh rule needs to be used to calculate correctly 
the velocity in the physical plane, (see Milne Thompson (1969), Clements (1973a, 
1973b, 1977) and Sarpakaya (1968). 

Another approach at body representation uses surface singularity sheets, most 
notably vortex sheets. Any potential solution <f> can be expressed, following Greens 
Theorems, as the sum of the influences of the boundary singularity sheets. The 
resulting integral equation is discretized into a set of algebraic equations and solved 
numerically for the set of panel strengths (sources, doublet or vortices) by enforcing 
the tangency boundary condition at the finite number of collocation points (Smith 
(1962), Smith and Hess (1966), Hunt (1978), Lewis (1981, 1987)). 

1.5 The Full- Vortex Cloud method 

A relatively recent development is the so-called ’full-vortex cloud’ modelling which 
attempts to simulate the flow through the a solution of Navier-Stokes equation by 
an operator- splitting method , first proposed by Chorin (1973,1978). It has been 
developed extensively by Smith and Stansby (1987), Lewis and Porthouse (1981) 
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and Porthouse(1983). 

In this method, the vorticity equation is divided into a convective and and 
diffusive part and two equations are solved sequentially. 

Thus, the Navier-Stokes equation 

dq 




V P , ..T72, 


+ z/V 2 q 


is first converted into the two dimensional vorticity transport equation 

du 


dt 


+ q.Vu; = i/V 2 u; 


and then is split into two equations 

du 

~di 


+ q.Vu; = 0 


and 


du 

dt 


= uV 2 u 


(1.3) 


(1.4) 


(1.5) 


( 1 . 6 ) 


or its non-dimensionalized version 


du 

~dt 



(1.7) 


The two equations are then solved sequentially. A practical manifestation starts 
with the potential flow about the body calculated with vortex sheet replacing the 
body surface. This vortex sheet is then discretized into a number of vortices. The 
vortices are then made to diffuse according to equation 1.6 by using some kind of 
random- walk technique. After the diffusive step, the vortices are then allowed to 
convect with the convection velocity determined by the influence of all the vorticity 
elements, the free stream and the new vortex sheet required to satisfy the tangency 
boundary condition. After convection, the tangency condition is satisfied again, and 
the new sheet is ready to be shed, diffused and convected. 
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1.6 Present work 

The present work is in essence an extension of the work started by Parag Kumar 
(1993) and continued by Verma (1994). Parag Kumar calculated the flow past an 
inclined flat plate where the body representation was obtained through conformal 
transformation, and the nascent vortices were shed at points whose location was 
obtain through Kutta condition and whose strength was based on the shear-layer 
velocity obtain as the average convection velocity of the first four vortices in the 
wake immediately next to the plate. 

Due to obvious limitation of this scheme in more complicated geometries Verma 
(1994) adopted a vortex panel scheme for body representation. The vortices were 
shed from fixed separation points and the angle of the separation streamlines was 
optimized through trial and error. Verma calculated the flow past a half cylindrical 
cup facing the flow representing a rather guide 2-D canopy. Verma also started, 
work on an operator-splitting technique for solving the Navier-Stokes equation and 
obtain some elementary solutions for flow past a cylinder. 

The present work takes off from the work of Verma and develops further the 
full-vortex cloud scheme involving operator-splitting. The first part of the work is 
concerned with exploration of the method as applied to the problem of boundary 
layer on a flat plate. This problem affords the opportunity to explore the problem 
of large Reynolds number flows and the use of distinct time steps for diffusion and 
convection. 

Next the full-vortex cloud analysis is developed for the problem of flow past a 
cylinder. The focus of the program was to develop, if possible, the capability to 
predict drag on bluff bodies with a view to evaluate changes in shapes intended for 
drag reduction. For this purpose, attempts have been made to develop a procedure 
for calculating the flow past multiple two-dimensional objects. 



Chapter 2 


SIMULATION OF A 
BOUNDARY LAYER 

2.1 Introduction 

As a preliminary to the development of full- vortex cloud method of calculating the 
viscous separating flow past a bluff body, let us first consider the computation of 
the development of a boundary layer on a flat plate by a discrete vortex method. 

The vorticity transport equation for the case of incompressible 2-dimens:onal 
flow is obtained by taking the curl of the Navier-Stokes equation to obtain 

^ + q ,Vw = i/V 2 u (2.1) 

where the vorticity intensification term cu.Vq has dropped out because of the two- 
dimensionality of the flow. 

In its non-dimensional form, this becomes 

q - v “ = -- 2) 
where Re stands for the Reynolds number. 
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It is convenient to consider the flow past a fiat plate (see Fig. 2.2) as a potential 
flow in which the no slip condition at the plate surface modelled as a source of 
vorticity there. 

The vorticity thus introduced into the flow then diffuses and convects according 
to the Eqn. 2.2 above. We will develop here the operator-splitting method which was 
first developed by Porthouse and Lewis(1981) in which the convection and diffusion 
instead of acting simultaneously are made to act serially, so that equation 2.2 is split 
into two equations 


and 


du 

dt 

dui 

~di 


-q .Vw 


Re 


(2.3) 

(2.4) 


Equation 2.3 is quite straight formed to handle, on reorganization, one obtain 


^ + q.Vw = ~ = 0 (2.5) 

dt Dt K ’ 

This is taken to imply that the discretized vortex elements in the flow are 
convected as material elements at the local flow velocities. 

Equation 2.4, which is the classical diffusion equation is solved using a random- 
walk method, We first develop an algorithm to solve this equation after the procedure 
suggested by Chorin(1973). 


2.2 Solution of the diffusion equation 


Batchelor(1970) obtained the well known solution of the problem of diffusion of a 
concentrated vortex of strength P which is located at the origin at time t — 0. The 


subsequent distribution of vorticity in time and space is obtained as 

r 


U) 




Avvt 


-£ 

e 


( 2 . 6 ) 
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which represents the radial distribution of vorticity. Lewis(1981) solved the same 
problem by discretizing the concentrated vortex into a large number N of vortex 
elements at time t = 0, all located at the origin, and then letting the various 
elemental vortices undergo a random walk. The strength of each of these elemental 
vortices is taken as T/N . 

If the solution after time t is given by Eqn 2.6, then the number n of elemental 
vortices in an area rdQdr at that time should be given by 

—— = uj(r,t)rdddr = rdQdr (2.7) 

N Aicvt v ' 

If p denotes the probability that a specific element ends up in this area element 
after time t, then 

n 1 =I ^_ 

p = — = e rdQdr (2.8) 

N 4xi A ' v ' 

An appropriate random-walk method then displaces each element i in the radial 
and angular direction by amounts and 6 t over the time interval 0 to t, such that 
the radial probability distribution given by Equation 2.8 is satisfied for all r_\#Ar 
elements which make up the given target area. 

The angular symmetry of the governing diffusion equation suggests that the 
scattering in the 6 direction should taken place such that an element has equal 
probability of ending up in any direction. This is achieved by generating a random- 
walk Qi between 0 to 1, and letting the element take up 

di = 2 %Qi (2.9) 

at the itfttime steps. 

To obtain the radial scattering of vortex elements first integrate Eqn 2.8 for 6 
between 0 and 2x and to obtain the probability that a given element will lie between 
r and r + dr, by 

1 =sL j 

n = e ivt rar 

2 id 


( 2 - 10 ) 
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This equation may next be integrated over dr for r between 0 and r to obtain 
the following equation 

P = 1-eTsr (2.11) 

where P gives the probability that the given vorticity element will lie somewhere 
within the circle o/ radius r. 

At this point it is interesting to notice that this probability agrees with the 
normal distribution curve of statistical theory. For the random- walk process, then 
we generate a random number P t between 0 to 1 with uniform distribution, Then 
we postulate that the ith vortex element will be found at a radius r t - such that 

-r? 

Pi = 1 - ( 2 . 12 ) 

The radial shift of each of the vorticity element can then be obtained if the above 
equation is inverted as follows 

n = y4 vtln(y) (2.13) 

2.2.1 Diffusion of a vortex 

Consider a vortex of strength 1.0 situated at the origin. To obtain the vortex 
diffusion with time, we implement the above scheme in the manner shown in the 
flow diagram which is shown on the next page, 

2.2.2 Random number generation (Step 1) 

Random numbers P { and Q { are generated by using formula. 

Pi = Pi + 1 — integer(exp(5.0ln(Pi + i + 1.10318))) (2-14) 
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by defining a series of annular bins separated by equally spaced the process, e.g. P 
= 0.5. The quality of random number has been verified by using statistical methods. 



2.2.3 Vortex diffusion over a time steps At (Step 2) 

Now consider the diffusion over a small time increment At. The displacements of 
element i during time At will then be 

Adi = 2tt Qi (2.15) 

A ri = y^ANn(-i-) (2.16) 

Thus after the increment At the new coordinate location as the i th element will 
become 

x'i = Xi + Ar; cos Adi 
y'i = yi + Ar,- sin Ad { 


(2.17) 

(2.18) 
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2.2.4 Radial distribution of vorticity 

The evaluation of vorticity for the scattered vortex elements, the process, e.g. P = 
0.5. radius rj. Suppose that rij elements are captured by the area lying between rj 
and rj + i . Thus, the vorticity at the jth strip will be estimated by 


cu(r) = 


n i 


(2.19) 


N vr(r? +1 - r]) 

For example, a unit strength vortex which was represented by 1000 elements of 
strength 0.001 and kinematic viscosity was 1.0, is given 10 random-walk over the 
time period t = 1.0 based upon the equation (2.15 to 2.18). The output including 

the comparison with the exact solution has been given in Table 2.1. 

4 jxy - - — 


200 - 


A 

I 

{ 000 — 

>. 


°o 


<*> 


'9 o a 

O o o 


- 200 — 


-4XO- 


x — — L.. — . ...i... J 

-4J0Q -2X30 0X30 2X30 4J0O 

X — > 


Fig 2.1. Random diffusion of 50 vortex elements over four time steps. 
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Table 2.1. Diffusion of a point vortex 


Bin no. 

No. of elements in bin 

rms rad. 

Numerical soln. 

Exact soln. 

1 

43 

0.282843 

0.085546 

0.078002 

2 

110 

0.632456 

0.072946 

0.072005 

3 

157 

1.019804 

0.062468 

0.061358 

4 

181 

1.414214 

0.051441 

0.048266 

5 

163 

1.811077 

0.036031 

0.035048 

6 

125 

2.209072 

0.022607 

0.023494 

7 

100 

2.607681 

0.015303 

0.014537 

8 

48 

3.006659 

0.006366 

0.008304 

9 

36 

3.405877 

0.004213 

0.004379 

10 

21 

3.805260 

0.002199 

0.002131 

11 

5 

4.204759 

0.000474 

0.000958 

12 

10 

4.604346 

0.000865 

0.000327 

13 

1 

5.003998 

0.000080 

0.000152 

14 

0 

5.403702 

0.000000 

0.000054 

15 

0 

5.803447 

0.000000 

0.000013 

16 

0 

6.203225 

0.000000 

0.000005 

17 

0 

6.603030 

0.000000 

0.000001 

18 

0 

7.002857 

0.000000 

0.000000 

19 

0 

7.402702 

0.000000 

0.000000 

20 

o 

7.802564 

0.000000 

0.000 000 


2.3 Discrete vortex simulation of the boundary 
layer 

When a fluid motiron is started impulsively from rest, vorticity is generated on the 
surfaces of all solid boundaries within the flow. From that moment on the viscosity 
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in the fluid gives rise to viscous diffusion of the vorticity in a direction normal to 
the surface of the body. If the surface vorticity is divided into discrete elements 
which are given random motions to simulate viscous diffusion and are subjected to 
convective influence of the free stream and all other vortices present in the flow, it 
would appear that a method is immediately available for calculating the development 
of a boundary layer flow. 

But the convection should be handled in a manner which models the presence 
of the flat plate as in the previous boundary correctly. This is handled here by 
introducing in the flow field mirror images about y = 0 of all the vortices present 
in the flow. The next effect of the presence of the image system with the system 
of the physical vortices is to make y = 0 a streamline (with zero normal velocity 
everywhere) and this is what is sufficient for maintaining the boundary condition. 
Fig. 2.2 shows an appropriate convection scheme based on the operator-splitting 
procedure discussed above 


U(x) 



Fig. 2.2. Mirror image system for modelling vorticity creation and convection in 

discrete vortex boundary layer simulation. 

The procedure is based on creating and shedding at the plate surface y — 0 a 
vortex sheet, of discretizing this vortex sheet into a number of line vortex elements, 
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and then to diffuse and convect. These vortices and other vortices HYe in the 
flow field appropriately. At each time step, the number of vortices are introduced at 
the plate surface. The development of the steady boundary-layer is then achieved 
through an unsteady process. 
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2.3.1 Vorticity creation and shedding (Step 1) 


According to the fundamental principle of the vortex-sheet method, the strength of 
the vortex sheet that needs to be postulated at a surface is equal to the slip velocity 
at that surface under the action of the inviscid potential flow and the vortices that 
have been shed up to that time. 

Thus, at time t = 0, the slip velocity is U all along the flat plate and the initial 
vortex strength is 7 (s t -) = U at y = 0 all along the wall surface 0 < x < l. Due 
to viscosity this sheet will immediately begin to diffuse in the manner discussed in 
section 1.2, and also undergo convection. At any later -L'm-e. the new vortex sheet 
at the plate will have strength given by 


N 


7 («*) = Ui = U - - 


AI >; 


* (*< - x i) 2 + yj 


( 2 . 20 ) 


where ATj is the strength of the discrete vortex j in the flow field (shed previously) 
and located currently at ( Xj,yj ). 

This vorticity is shed freely from the surface and is represented by M new discrete 
vortices of strength AT,- = 7 ( 3 , -)As,', now free to undergo diffusion and convection. 

In calculating the slip velocity, the summation in Equation 2.19 has to be carried 
out over all vortices and their mirror images, which are introduced to take care of 
the normal velocity boundary condition. 

2.3.2 Viscous diffusion (Step 3) 

After the creation and shedding of the new vorticity, all other discrete vortices must 
be subjected to random walks to simulate viscous diffusion during the time step 
At. A complication that results in this procedure that about half the newly created 
vortices, and some of the other vortices, diffuse to the half region y < 0. This clearly 
is unacceptable because of the impervious wall boundary condition. Chorin (1978) 
used a scheme in which he started with a double strength 27 ( 3 ,- )As^ with about 
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one half of the vortices that are diffused under the plate being simply quenched. 
This results m about the right amount of vorticity being introduced in the flow. 

Porthouse (1983) recommended that the nascent vortices may be introduced 
at an offset from the wall. Based on the consideration of the root mean squared 
development of vorticity, he recommended a displacement e = -jAi/At/3. 

We have used the Porthouse recommendation and in addition used a reflection 
procedure interior of the body (below the plate y = 0, 0 < x < 1 in this present case 
of the flat plate) is reflected back normally the same distance. This conserves the 
vorticity in the diffusion as well as the convection operation. 

2.3.3 Vortex convection (Step 5) 

As we have noticed the governing equation for the vortex convection phenomenon 
is 

du 

— + q.Vui = 0 (2.21) 

A numerical simulation of this equation suitable for computation is possible if a 
large number of point vortices are taken to simulate the area distribution of vorticity 
within a given region of the flow. Each element will then move under the i nfl uence 
of all the other elements. Vortex convection in this simple geometry of flat plate is 
correctly modelled by the mirror image system as shown in Fig 2.2. The convection 
velocity components experienced by element ATj due to element AT; and its mirror 
image are then given by 

A Uj = 

A Vj = 

where, 

r l = ( X J “ X if + (&' _ Vi ) 2 



(2.24) 
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r 2 ~ ( x j ~ x i ) 2 + (j Ij + Vi) 2 


( 2 . 25 ) 

( 2 . 26 ) 


Now, we can sum up the influence of all the vortices in the cloud as 

i=N 

Uj'doud — 'y 1 

i=l 

*W 

i=N AT- / 1 1 V 

Vj, cloud = 2 tt ~ Xi ^ \ r% ~ ^2 J ( 2 . 28 ) 

*¥i 

In addition to this we have the self-induced velocity of element AT; obtained through 
the effect of its mirror image vortex. 


Ar » ( Vj ~ Vi _ yj+JJj 
2 tt \ r? To 


( 2 . 27 ) 


Ar- 

Au,- = , Aui = 0 (2.29) 

The resultant velocity of the i th vortex is simply given by the combined influence of 
all the terms and the external velocity. 

— Ui, external Ui, cloud Au,' (2.30) 

= Ft, external "4" cloud (2.31) 

Applying a forward difference for the elementary time step At the new position of 
the vortex is then computed from 

WSU* ■ . 

x\ = + Ui^t (2.32) 

■y'i = Vi + u «At (2.33) 

Since all irreversibilities arise because of viscous diffusion which is not included in 
this step, the convective algorithm is developed so that it is completely invisible. 
This has been ensured by using a multi-step prediction- correction scheme. 
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2.3.4 Merge of vortices in close proximity (Step 6) 

Th ough the discrete vortex modtfl gives fairly good representation of the flow fields, 
the nature of modelling in which region of viscous flow with vorticity are replaced 
by inviscid free vortices, with velocities tending to infinite as one approaches the 
line vortex, makes the flow field close to a vortex to be highly vortex. This results 
in the convection velocities of vortices close to other vortices to be insignificant and 
random errors. 

Various schemes have been tried to overcome this problem, and all of them 
depends upon smudging a point vortex into a vortex blob, and any blobs which 
touch one-another are merged algebraically. 

We have developed a procedure in which each newly created sheet vortex element 
of strength 7(3.) and element length As, is immediately redistributed as a Rankine 
vortex of radius r 0 = A s n / tc for which the induced velocity at radius r is given by 

ATr 

Ag = 2^2 f or r < r ° (2.34) 

= 7 T~ for r>r 0 (2.35) 

2nr 

Any vortices which drift close enough for their vortex cores to be overlapping are 
then merged into a single vortex and the resultant strength of the vortices will be 
algebraic sum of the individual vortices and it will be placed at the center of vortices. 
For example two vortices of strength I\ and T 2 located at positions (xl,yl), (x2,y2) 
which happen to be closer to one another than a specified limit r min , which is equal 
to the sum of the radius of the two interacting Rankine vortex. 

Thus, vortices which axe drifted closer than the r min are merged into a single 
vortex and after merging the two vortices the strength of the resulting vortex is 
given by 

(2.36) 


r = r*i + r 2 
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at the position 


l£ij x i + jr 2 [ a: 2 

" l r il + |r 2 | 

l£j_l_yi + 1^2] y 2 

l^il + |r 2 | 


(2.37) 

(2.38) 


The combination produces two results. First it ensures that unrealistically large 
convection velocities are not introduced, and second, that the total number of 
vortices in the field remains manageable. 


2.3.5 Calculation of velocity profile (Step 9) 

The velocity profile is calculated by integrated the vorticity profile, neglecting the 
velocity component v in the y direction, which is a standard boundary -layer theory 
assumption. Thus, 

u (y) = / u(y) d v (2-39) 

Jo 

Numerically, the co(y) profile is obtain by establishing at the required location x, 
sampling boxes of width Ax and height Ay spacing the total number of vorticity 
elements present in each boxes at a given time step are counted, and then u.-(y,) is 
obtained by dividing the box at height y,- by its area AxA y 

Then, the velocity at y n is obtained by 

0 


ryn 

(yn) = / u(y)dy 

Jo 

(2.40) 

n 

(2.41) 


i-1 


2.4 Results and Discussions 

A typical vortex cloud, due to the scatter of vorticity elements, as it develops on the 
flat plate is shown in Fig 2.3. 60 surface elements were used to represents a plate 
of length l = 1.0 with a uniform mainstream flow at U = 1.0 Fig. 2.4. shows the 
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velocity profile as obtained by the above scheme for Re 500, 1000 and 2000. The 
Blasius profile is also shown in the figures for comparison. 



Fig 2.3. Laminar boundary layer on a flat plate with constant mainstream 
velocity U 


Selection of element size and time step 

An appropriate device of the time step A t is one in which the convective and diff usive 
displacement are of the order of the element size. Thus, the convective displacement 
6 C which is of order |£/A£ C should be of the same order as the element size, say 
M being the number of panels taken along the plate. This gives the convective time 
step as A t c « 

The typically diffusive displacement on the other hand is 

Ar,- = [4i/ A f M-f)] 1/2 (2.42) 

obtain from Equation 2.25. by as 8 d rZj [kvAtln* 2] 1 / 2 , so that 

8 d » [4:uUn2/MU} 1/2 (2.43) 
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or 



(2.44) 

(2.45) 


or 


At = 


16Un2 
U Re 


(2.46) 


For very large Reynolds numbers then, we see that the diffuse displacement time 
step is extremely small. A considerable saving in calculation effort can be effected 
by choosing diffbrcwt time scales for diffusion and convection. 

Thus, with reference to Equation 2.14 , we choose At = Atp such that 
so that 


A to — 


P 

l6uM 2 ln2 


(2.47) 


The ratio of the diffusive time step Afp to the convective At c then is 


A t D _ _ Re 

At c * 16M7n2 


(2.48) 


For a value of Re= 10 s and M — 100, we obtain a value of N t = 90. This means 
that for each 90 convective steps we need for perform only one diffusive step. This 
results is a very considerable saving in time of computation. 

Fig. 2.5 shows the boundary layer profile generation by vortex dynamics for Re 
50000 and 100000. At these Reynolds number, most boundary layer are still laminar, 
therefore close to turbulent. But the calculated profile is quite different from the 
Blasius profile suggested that the numerical turbulence implied by the finite scale 
As / t r of the discrete vortices perturbs the shear layer leading to transition. 
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Vortices close to walls 


The vortices which are too close to wall, develop the boundary-layer velocitv profile 
according to the law of the wall, 


u 

u T 


yu r 

V 

2.44/n 



+ 4.9 


for** < 11 


for** > 11 
v 


(2.49) 

(2.50) 


where the friction velocity u r is defined in term of wall shear stress tq by 



(2.51) 

(2.52) 


Thus u T may be directly related to the vorticity w(y 0 ) close to the wall. The suitable 
approach here to calculate a>(y 0 ) is to take the average vorticity strength over the 
inner half of the sub-layer. The logarithmic scale of plotting then helps to draw out 
the characteristics of the boundary layer in both the laminar near wall and turbulent 
outer regions. The solution for Re (50000 and 100000) has been shown in Fig. 2.6. 
The results which comes directly from the law of the wall are also incorporated for 
the comparison. 




— 

Blasius solution 


c 

o 

Vortex sol. for 

Re=500 

□ 

D 

Vortex sol. for 

Re=±000 

A 

A 

Vortex sol. for 

Re=2000 


Fig, 2.4. Laminar boundary layer on a flat plate at Re= 500, 1000 and 2000 












Chapter 3 


THE FULL- VORTEX CLOUD 
ANALYSIS OF FLOW PAST 
BLUFF BODIES 

3.1 General Introduction 

In this chapter the method of the full-vortex cloud developed in the last chapter is 

extended for flow past two-dimensional bodies of arbitrary shapes. The method is 

c 

an extension of the traditional discrete vortex methods which are valid for infinitely 
large Reynolds number (th ough with the postulation of the separated wake), and 
attempt to include the effect of the finite (but large) Reynolds numbers, and is based 
essentially on the procedure developed by Verma (1994). 

The surface of the two-dimensional body is covered with discretized vortex- 
panels. The strength of which is determined using Martensen method which is 
based on satisfying the zero tangential velocity condition by postulating a vortex 
sheet at the body surface. 

Following the treatment developed in the last chapter, this vortex sheet is seen 
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as the incipient boundary layer, which develops as this vorticity diffuses into the 
flow field and is convected downstream. The numerical scheme A developed in this 
chapter. Th^ basic strategy is displayed in the given flow diagram and following the 
operator splitting procedure developed in the previous chapter. 


3.2 Martensen method 


Consider the flow past a two-dimensional body of an arbitrary shape in the (x-y) 
plane, immersed in a uniform flow at an angle of attack measured with 
respect to x-axis as shown in Fig. 3.1. We may discretize the body surface into M 
elements. Let the length of n th panel be ds n and let it be covered bv a vorticity 
sheet of strength 7 (s) per unit length. The distance s is measured in the clockwise 
direction around the body perimeter from some datum which is usually chosen as 
the leading edge of the body. 


v> 



Fig. 3.1. Discrete surface vorticity modal for a 2-D body 


The velocity dq {j induced at s t due to a small vorticity element 7 {sj)dsj located 
at Sj elsewhere on the body directly follows from the Biot-Savart law, which in this 
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7 {sj)dsj 
27 T ry 


case reduces to: 

% = ~ ^ - (3.1) 

where r^- is the distance between the elements i and j. Now resolve dq^j parallel to 
the body surface at i where the profile slope is expressed as We express dq^ in 
terms of local coordinates by introduce a figure here to show coordinates first. 




(3.2) 


dVa 


-7 {sj)dsj 


-(*» - x j ) 


Resolving dqij parallel to the surface element s,-, we obtain the following tangential 
velocity component 

dvij = dUij cos /?,- + dVij sin /?,- (3.4) 

Integrating over the entire surface and including the effect of the non-singular 
Fredholm integral equation at the element $,• 


0 . 57 ( 5 ^) + k(si,Sj)j(sj)dsj + Woo (cos »oo cos fa + sinaoo sin /?,-) = 0 (3.5) 


Here the contour integral yields velocity parallel to the i th element due to all 
bound vortices around the profile. 

In the above equation, the last term is the component of resolved parallel 
to the body surface at i and k($i,Sj) is the coupling coefficient given by 


k(. s ii s j) 


1 ivi - yj) cos $i - (Xj - Xj) sin ft 


(3.6) 


2 tt (xi - Xj) 2 + (yi - yq) 2 ■ 

For convenience we may write a modified coupling coefficient K (s,-, Sj) by ab- 
sorbing 0.57(s,-) into it. Thus, i.e. 


K(si, sj) = k(si,sj)dsj - 0,50 Sijdsi 


(3.7) 
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Hence equation (3.5) may be written in the numerical form as 

M 

J2 K ( 5 i> s jM s j) = -(Uoo COS A + Ko Sin A) (3.8) 

7 = 1 

This is the basic equation of Martensen procedure. The$e are the details of the 
numerical strategies adopted here. 


3.2.1 Surface element geometry 

A set of (M+l) input data coordinates ( x n ,y n ) are specified as shown in Fig. 3.2. 
moving clockwise around the profile from the leading edge. Point M + l coincides 
with point 1, ensuring profile closure. Straight panels are obtained by joining 
successive data points. Elements lengths are then given by 

As n = ^/(x n+1 - ar n ) 2 + (y n+1 - y n ) 2 (3.9) 

Profile slopes and pivotal points which are located at the center of each element, 
are then given by 


COS fin — (^n+1 ^n)/ASj 

(3.10) 

sin Pn = G/n+l - J/n)/Asj 

(3.11) 

2-n = 2 (^"n+1 ”b ^n) 

(3.12) 

Vn = yiVn+l + Vn) 

(3.13) 



33 



Fig3.2 Surface element geometry 


3.2.2 Determination of Coupling coefficients 

By definition, coupling coefficients K(s,-, Sj) represents the velocity at Si parallel to 
the body surface induced by a unit vortex located at the j th panel. The determina- 
tion is rather straight forward except for some complications as before. 

Effect of Curvature 

A curved panel has non-zero value of the self-induced velocity as shown in Fig. 
3.3. But the discretization scheme adopted uses straight panels which necessarily 
have zero self-induction. For better accuracy, the coupling coefficients must be 
modified so as to take into account the curvature of the body surface. Generally, 
such calculation involves the use of curve fitting methods, which are not suitable for 
the cases in which the curvature varies enormously and rapidly in certain locations 
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of the body profiles. Lewis gives a simple and alternate relation which has been 
used in the programs here, as 

A'(s,-,»,-)=0.5- ^'+ 1 8 ~ ft -' ) (3.14) 



(b) Finite seif-induced parallel velocity 
of a curved vorticity element 

Fig. 3.3. Self-induced velocity of a surface vorticity element. 

Use of sub-elements 

When the body is thin, the elements opposite a given collocation point come so 
close that the concentration of the distributed vorticity at the middle point may 
result in overestimating the influence of that element. This result is manifested 
mathematically, by the growth of the back diagonal terms in the matrix of coupling 
coefficient to the order of the diagonal elements (approx 0.5 ). The easiest way 
to overcome this is to breakup such elements into a number of sub-elements and 
to replace the relevant coefficient by the average of the coupling coefficients of the 
sub-elements. ' 

Consider a case in which a element Sj is broken up into N subelements each of 
length dsj/N. The coupling coefficient of this the i th element then is given by 

K(s h y.Jj- K( Si , s^dsj (3.15) 
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or in the numerical form 


where 




I N 

N 


k=l 


(y» ~ Vk) cos /3j - (z,- — x k ) sin /?,- 
( x i ~ x k ) 2 + {Vi - Vk ) 2 


*k = Xj + {k- 0.5 (N + 1 

y t = y. + (k- 0.5 (N + l)) Aj, ' Smft 


(3.16) 


(3.17) 

(3.18) 


It is found that use of two or three sub elements gives better results than using 
more subelements. We shall use only two subelements in all our calculations. 


3.3 Martensen analysis in the presence of a cloud 
of vortices 

In a quasi-steady flow the potential flow past a 2-D body of arbitrary shape in the 
presence of a cloud of line vortices is described by the boundary integral equation 


-0.57(s m ) + )l(s n )ds n + (Uoo cos J3 m + Voo sin #) 

z 

+ ATj(U m j cos /5 m + V mj sin (3 m ) = 0 
i= i 


(3.19) 


where unit velocities U m j and V m j are given by 


U, 


1 {Vm ~ Vs) 


m3 


5 Knj ~ * 


1 (x m Xj ) 


2tt' rlj — . 27T r* 


(3.20) 


with 


T mj — ( Xm X j) (^ m y <?) 


Here the last term representsthe contribution to the Dirichlet boundary condition 
of zero tangential velocity of the Z line vortex elements of strength Fj. In the scheme 
of calculation developed here, the last two term are shifted to the RHS y at each 
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step the matrix of influence coefficients does not change while the RHS changes (due 
to convection of vortices and introduction of nascent vortices), it is convenient to 
use the matrix inversion path to solve the equation. 

3.4 Vortex shedding from the body surface 

In the application of discrete vortex modelling for the separated flow about bluff 
bodies, the new vortices are shed into flow field at each time step. Just as in 
the development of boundary layer described in last chapter, the surface vorticity 
distribution required to satisfy the no-slip condition at the body at any time step is 
taken as the vorticity created in the time interval, and as the vorticity diffuses into 
the flow field. And as before, this vorticity diffuses and convects (essentially along 
the body surface) forms the thin boundary layer and the wake of the body. 

As was asigned in the last chapter, if the diffusion of the vorticity is taken to 
start from the location of the vortex surface, i.e. at the body surface itself, about 
one-half of the vortices will move into the body interior right at the first step. Chorin 
suggests using a vortex sheet of twice the strength so that after the first diffusion, 
we get the correct strength of the vortices. 

An alternate approach is the release of vorticity at some offset from the surface. 
The offset is determined by equating it to the root mean square drift of vortices in a 
random-walk step, which is yfevEt/Z. In fact, this step is identical in result to the 
Chorin scheme, except that the diffusive calculation are reduced by a factor of 2. 

3.5 Diffusion of vortices 

Diffusion of vortices is modelled by a random walk method similar to that used in 
the last chapter. Here unlike Chorin(1973), Lewis(1981) and Verma(1994), we have 
used the point vortex diffusion model, walking randomly in the 6 direction with flat 
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probabily, and walking randomly in r direction with normal distribution, as shown 
in the last chapter. The previous researchershave permitted diffusion only along the 
body surface on the plea that the resulting term in the boundary layer equations 
axe l^jr. But if we note that the original term was V 2 cu from with the x direction 
was neglected(compared to x convection), it is clear that using the simpler point 
diffusion model is quite correct with the added advantage of the formulation already 
available. 


3.6 Convection of vortices 

Once the vortices are shed from the body surface, they are free to undergo convection 
under the combined influence of vortex cloud, bound vorticity on the body surface, 
and external flow field. 

Let us focus out attention on a single vortex ’m’ located in the flow field at the 
position (x m ,y m ) Fig. 3.2. If there are Z vortices in the cloud, the velocity at m th 
vortex will be influence by these vortices. If we consider n th vortex in the cloud 
located at the position (*„, f /*), the velocity at m due to the vortex n will be given 
by Biot-Savart law as 

rr F n (y m ~ Vn) ( Q 91 >1 

U mn = — • o 


Bn (x„i x n ) 
27r' rl 


we can sum up the influence of Z vortices as 


U mi doud 


Tn f IJm Vn 


( 3 . 22 ) 


( 3 . 23 ) 


r„ f x m Xn 

V m , cloud 2—J | „2 


( 3 . 24 ) 
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If there are M discrete panels on the body surface with distributed vorticity of the 
strength 7 (s n )ds n at the n t/l element, the combined influence of these bound vortices 
on the m th free vortex in the wake will be given by 

'y(sj')dsj ( y m Vj\ 


U, 


m, bound 


E 

3=1 


2tt 


mj 


V _ _ 7 i s j)dSj f x j 

Vm, bound- 2 7T \ r 2 mj 


(3.25) 


(3.26) 


The velocity far upstream the body influence the velocities of the free vortices as 


Um, external ^oo 

(3.27) 

Vm, external ^oo 

(3.28) 


The resultant velocity of the m th vortex is simply given by the combined influence 
of all three terms. Thus we can write 


Uyn — external T U m , cloud ~b u m ,1 bound 


(3.29) 


V m — Trn , external T Kn,i cloud "b ^m, bound (3.30) 

Now, for applying the vortex convection, Let us consider a vortices at the position 
of a and then by using forward difference approach, the revised locations of the free 
vortices during the time interval is given by 


%mb — %ma 

' 

(3.31) 

Vmb == Vma + 


(3.32) 


After convectiwj all vortices, if we now recalculate the new convection velocities at 
locations b, (u mb , v mb ), we may march forward again through next time step to the 
final position. For better estimation of vortex convection, the central difference 
technique in which, the effect of the curvature of the drift path is also considerable 

is used like that, 
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•Em, final Ima H~ 0.5(u ma -f U m j)At (3.33) 

Vm, final = Vma + 0.5(u ma + V m b) At (3.34) 

3.6.1 Convection of vortices in very close proximity to 
body surface 

In the last section, we had mentioned some difficulties while calculating the coupling 
coefficient by the use of Biot-Savart law alone. Such problems are again encountered 
when calculating the velocities of vortices which have drifted too close to the body 
surface. It is observed that even the use of sub- elements does not yield good results 
when the gap ratio e/A s is less than 0.40. In such cases we go for the ’Mirror Image 
Technique’ to calculate the velocities. 


Mirror image technique 


If the gap ratio e/ As is less than 0.4Q^ then the velocity of the vortex is calculated by 
the use of simple mirror image technique. If we assume that the panels are straight, 
then the convective velocity of the vortex is given by 


<ld = 


ATn 

47TC 


(3.35) 


where qd is parallel to the panel. While using this formula we neglect the Uoo and 
Ko components of the velocities. 

We shall summarise the techniques available to us for the calculation of velocities 
of free vortices 


(1) For gap ratios, e/As > 1.0 normal surface vorticity modelling will suffice 
without the use of sub-elements. 

(2) For 1.0 > e/A s > 0.4 sub-elements should be used when calculating the 
convective velocity due to the bound vorticity. 
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(3) For e/ As < 0.40 mirror image modelling should be used. 



Fig 3.4. Mirror Image system for self-induced convection 

Some other difficulties that arise when simulating |lon> behind bluff bodies using 
vortex element techniques are 

1 Vortices may accidentally stray inside body contour due to inaccurate convec- 
tion routines. 

2 Vortices may drift too close to the body (so that even mirror image technique 
fails to calculate the actual velocity). These vortices may generate s«ricms potential 
flow errors or excessive self-convection velocities. 

Selection of element size and time step 

An appropriate device of the time step At is one in which the convective and diffusive 
displacement are of the order of the element size. Thus, the convective displacement 
S c which is of order At c should be of the same order as the element size, say , 

M being the number of panels taken along the body. This gives the convective time 
step as A t c « jjjj. 
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The typically diffusive displacement on the other hand is 

Ar t - = [4i/A (3.36) 

•*» 

got from Equation 2.25. by as 6 d w [4^At/n2] 1/2 , so that 

6 d « [4idfn2/M[/’] 1 / 2 (3.37) 


or 


2 M KUl/ 

(3.38) 

n 6My 

V Re J 

• (3.39) 

16//n2 

At = 

U Re 

(3.40) 


For very large Reynolds numbers then, we see that the diffuse displacement time 
step is extremely small. A considerable saving in calculation effort can be effected 
by choosing difficrcwt time scales for diffusion and convection. 

Thus, with reference to Equation 2.14 , we choose At = Atp such that j, 


so that 

A ^ 

At ° " 16vM 2 ln2 

The ratio of the diffusive time step A t D to the convective A t e then is 

—2. _ at = Re 
At c “ * !QMln2 


(3.41) 


(3.42) 


3.7 Zero circulation condition for RHS term 

The necessary condition that the net circulation around the cylinder perimeter due 
to the externally located vortex T must be zero, requires that, all the vortices in the 
wake behind a body should produce net zero circulation on the body surface. If this 
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condition were not satisfied, there will be an implied erroneous residual vorticity 
bound within the body profile leading to a flux through the body surface for this 
shown a zero circulation condition is built into the RHS. The influence of the vortex 
cloud on the m element on the body surface makes that following contribution to 
the RHS 

. .. z 

RHS miC i oud = £ ATj ( U m j cos 0 m + V m j sin /5 m ) (3.43) 

j= 1 

Circulation around the whole body surface can. then be written as 

^ / RH Sn^cioudASn ^ 0.0 (3.44) 

n=l 

To ensure this, then, at each step the largest value of RHS ? , c i oud is replace by 

M 

RHS PiC i oud = — — — Y2 RH Sp,ci 0 ud&Sn (3.45) 

^ S P n=l 

n/neqp 

3.7.1 Conservation of vortices 

In accordance with Helmholtz theorem the vorticity should be conserved at all times 
and hence if the vorticity is zero everywhere when the flow starts, it should be zero 
all through. We can write this in terms of mathematical formula 

M z 

£ 7 {Sn)ds n + } y AT j Voire — 0 (3.46) 

71=1 j- 1 

where r ci v c is the cumulative strength of all the vortices which are snuffed out for 
various reasons as the calculations proceeds. This equation must be added to each 
of the governing equations to ensure that vorticity is continually being conserved. 
Thus, we can write the most general equation as follows. 

M 

(k[s mi s n ) + As n )7 («s n ) = —Woo. cos(o: 00 — (3 m ) 

7X=1 

z 

+r circ - £Aiy(i + i7 mi COS /5m V m j sin /5m) 

j=l 


(3.47) 
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3.8 Merge and Reflection of the Vortices 

As we are releasing M vortices at every time step the number of vortices in the 
flow field will grow unless some action is taken. There is a need to curb the rise in 
number of vortices in the flow field because of the following reasons 

1. The CPU time is directly proportional to the N 2 , where N is the number of 
vortices in the flow field. 

2. As the vortices come closer, they will induce excessive velocities on each other. 

3. The influence of the vortices on the body surface decreases with the distance. 
Therefore computational efforts can be minimized if we can reduce the number of 
vortices in the flow field. 

For these reasons, we shall merge those vortices which have come too close to 
each other than some specified distance. The merging is done using exactly the same 
manner as done in the previous chapter for calculating the boundary layer on a flat 
plate. 

3.9 Calculation of surface pressure distribution 
and forces 

Once convection and diffusion have been completed in the numerical simulation The 
Navier- Stokes equations reduces to 

^ = -i Vp, (3.48) 

dtp 

At any point s n on the body surface, the velocity q parallel to the surface is 
given by 7 (s„). Therefore, 

d]Ps___ d7(sn) 


(W) 
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from which we can directly calculate the numerical for the change in surface 
pressure over the surface element n during the time step At value of 


Ap n = -p 


7pn)As n 

At 



( 3 . 50 ) 


Equation may be integrated to yield the surface pressure at any point on the 
body 


P M 

Pm=Psta g ---T T i ( 3 - 51 ) 

^ 1 J+l 

Average lift and coefficients can be calculated from the values of pressures. 
Spectral analysis technique is used for finding out the dominant frequencies f of 
the excitation forces. From this dimensionless frequency can be expressed from the 
Strouhal number 



J± 

Woo 


( 3 . 52 ) 


3.10 Determination of Lift and Drag coefficient 


Once the static pressure distribution is known around the body surface, lift and 
drag coefficient can be estimated from the following equations 


M 


C L 


0.5 pWll pWl 


Y,Pn cos/3 n As n 


C D = 


D 


M 


^2 p n sin Asj, 


( 3 . 53 ) 

( 3 . 54 ) 


0.5p^/ pWl 

where 1 is the characteristic length and p n = 0.5(p n + p n +i) i s ^he average pressure 


on element n 
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3.11 Results and discussions 

For the application of full vortex cloud analysis we have considered two geometries, 
a circular cylinder and a square cylinder. We have studied the flow past these bodies 
at various Reynolds number. In these applications, we first divided the body surface 
into a number panel. Vortices of appropriate strengths are shed from each el em ents 
normal to the surface at each time steps. These newly created vortices axe subjected 
to diffuse by the random walk technique as discussed earlier in detail. After that we 
merge the two vortices into one vortex when the distance between these two vortices 
is less than the specified distance and vortices which drifted under the body surface 
are reflected back to that position at each iteration. The time step is taken around 
0.12 as per rule and the case was run for 250 time steps ( T max = 30.0) which is 
sufficiently large to cover several oscillations of the vortex sheet. Finally, the drag 
and lift coefficients has been calculated for both the bodies at various Reynolds 
number ( 20000, 30000, 40000 and 50000 ) and then, strouhal number has been also 
calculated by using the technique of spectral analysis. 

Flow past the circular cylinder 

For the case of circular cylinder in which the main stream flow is 1.0, we take 
the diameter of the circular cylinder 1.0 and discretize the body into the 32 surface 
element. For eliminating the numerical noise, we are averaging over two time steps in 
the calculation of lift and drag coefficients. The flow diagram for different Reynolds 
number are shown in Fig. 3.5. to Fig. 3.8. in which we have considered three 
different times ( 6.0, 12.0 and 30.0 ). where as the variation of the lift and drag 
coefficients w.r to time is shown from Fig. 3.9. to Fig. 3.13. Table 3.1. shows the, 
average values of the drag and lift coefficients with the Strouhal number. But there 
was lot of numerical noise in the pressure calculation so that we could get quite 
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higher values of drag and St. 

Table 3.1. For the case of circular cylinder 


Re 

Average 

Drag coefficients 

Average 

Lift coefficients 

St- 

20000 

3.1144 

6.655E-02 

0.3448 

30000 

3.2547 

4.384E-02 

0.3540 

40000 

3.2466 

1.424E-02 

0.3478 

50000 

2.8849 

-1.674E-02 

0.3288 


Flow past a square cylinder 

For the case of a square cylinder, main stream flow is same as circular cylinder and 
the side of the square is also 1.0 we discretize the body into 36 surface element. The 
flow diagram for different Reynolds number are shown in Fig. 3.14. to Fig. 3.17. 
as the case of circular cylinder and the variation of the lift and drag coefficients w.r 
to time are shown from Fig 3.18 to Fig 3.21. Table 3.2 is in same pattern as Table 
3.1, for the square cylinder. We get higher values of drag due to the problem of 

numerical noise in the pressure calculation. This is very difficult problem which we 

• \ 

face in this method. 

Table 3.2. For the case of square cylinder 


Re 

Average 

Drag coefficients 

Average 

Lift coefficients 

St 

20000 

3.6827 

3.324E-02 

0.2663 

30000 

3.6743 

4.526 E-02 

0.2690 

40000 

3.8223 

3.628E-02 

0.2681 

50000 

3.8013 

4.228E-02 

0.2672 
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( A) Tmax- : 6.0 



(B) Tmax~12.0 



( C) Tmax=-30.0 



Fig. 3.5. Full vortex cloud solution for circular cylinder at Re = 20000 



(A) Tmax--6.0 







(B) TmaX“12.0 




(C) Trnax=30.0 

Fig. 3.6. Full vortex cloud solution for circular cylinder at Re 


30000 



49 



(A) T max- 6.0 


CST* 


.Mr-;; * ^ 


* *.** , 
*** 

, . »* / * „ 


K w‘/*v* 


( 6 ) Tmax- 12.0 




‘ r* .u * x tJ * .. ‘Jfer- 


(c) Tmax-- 30.0 


c>% jft; \V # • "V' \:t \ 

' ’ - 0&' 

i- » * 


. * * «► 
fc « ir < j ! > % v* * < 


l ig. 3.7. Full vortex cloud solution for circular cylinder at Re = 40000 
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(A) Tmax-6.0 



fB') Tmax«i2.0 



:'C) Tmax-30.0 


Fig, 3.8. Full vortex cloud solution for circular cylinder at Re = 50000 
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(C) Tmax~30.0 



Fig. 3.13. Full vortex cloud solution for square cross section at Re = 20000 
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Fig. 3.14. Full vortex cloud solution for square cross section at Re = 30000 




( A) Tmax-6.0 



(8) Tmax-12.0 



(C) Tmax : -=30.0 


Fig. 3.15. Full vortex cloud solution for square cross section at Re = 40000 
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(A) T max=6.0 



(B) Tmax=±2.0 
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Fig. 3.16. Full vortex cloud solution for square cross section at Re — 50000 
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Jig. 3,17. Lift and drag coefficients on a square cross section at Re = 20000 j 
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Fig. 3.18. Lift and drag coefficients on a square cross section at Re = 30000 













Chapter 4 


THE FLOW PAST THE 
MULTIPLE BODIES 


4.1 General Introduction 

This chapter describes an attempt to simulate the flow past multiple two dimensional 
bodies by full- vortex cloud method. The procedure closely follows -Bitfc described in 
the previous chapter for a single body. 

Consider the flow past the P bodies of an arbitrary shape in the (x-y) plane, 
immersed in a uniform flow at an angle of attack measured with respect 
to x-axis. We discretize the surface qth into M q elements with q taking values from 
1 to P, Then, the flow past an assembly of P mutually interacting bodies may be 
represented by a simple adaptation of equation (3.8) to read 

p M q 

£ £ iC.(V ) = -(«» cos sin /?,„) (4.1) 

g= 1 7i=l 

where 

p = 1,2....P and m = 1,2 M p 
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Equation 4.1 states the Dirichlet boundary condition at surface element i of 
body p. In the case of multiple bodies the coupling coefficients matrix includes 
the contributions from all the elements ( n = 1,2 ...M q ) on each of the bodies 
( q = 1,2 ...P ). The number of elements M q chosen for each body may differ 
depending upon the individual geometrical requirements. 

The coupling coefficient representing the induced velocity at the pivotal point m 
of body p due to element n of body q is then given by 

jrpq f (|/pm Vqn) cos /^pm (^pm s i n Ppm 

1 mn ~ 2 ?r \ ( X pm - x qn ) 2 + ( y pm -y qn ) 2 

For this case the self - inducing coupling coefficient when (p — q and m = n) are 
given as before single body in equation 2.12. 



Kvv = n 6 - Priam jrimll (4.3) 

Jy mm v J 

The present investigation was started with a shape of finiting, the relative reduc- 
tion of drag by the various edge treatments that are being proposed for reducing drag 
of trucks and like shaped vehicles. The geometry for this example is a preliminary 
shape chosen to develop confidence in the method developed. We choose a square in 
a close proximity of a rectangle as shown in Fig. 4.1. Then, the general equation to 
solve the potential flow can be expressed in matrix form as follows where the right 
hand side remains same as discussed earlier. 

An Al2 1 (t( s )) = ( rhs ) ( 4 - 4 ) 

^ A2I A22 J 

where the coupling coefficient matrix has been partitioned to clarify the various 
body interaction. Thus typically, the sub- matrix A u contain all of the coefficients 
accounting for the interference experienced by first body due to the second body. 
The submatrix forming the main diagonal A n ,A 22 etc, account for the influence 
of each body upon itself and are thus identical to those obtained for each body 
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considered in isolation. It may be noticed that there are a large number of cross- 
influence terms signifying that the interactive effect is critical. 

The most critical aspect of formulation was the treatment of the terms involve 
panels which are in close proximity but happen to be on diffused bodies, and in 
keeping track of zero circular condition, and in effecting internal circular correction. 

4.2 Modification for bodies in close proximity 

When the bodies are very close to each other, there is the certainty of numerical 
errors should the gap between point pm and qn on adjacent bodies be less than the 
local element lengths A s pm or As qn , In that kind of case modification is needed to 
the i t h coupling coefficients in column n by the value 



m^i 

The k th element of body p should be the one in close proximity to element n of body 
q, as shown in Fig. 4.1. 


Fig. 4.1. Flow past a square in close proximity of a rectangle. 
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4.3 Results and discussions 

For the application of flow past the multiple two dimensional bodies, we have 
considered two geometries (a square cylinder and a rectangle cylinder) and, then 
studied at various Reynolds number (20000, 30000, 40000 and 50000). 

The body of the square cylinder and the rectangle cylinder were discretized into 
36 surface elements and 22 elements. The side of the square cylinder was taken 1.0 
and the sides of the rectangle cylinder were taken 1.0 and 0.2. Thus, 60 vortices of 
appropriate strengths are shedded from the bodies surfaces at each time step, and, 
then diffusion and convection were as usual. 

The flow diagrams are shown in Fig. 4.2 to Fig. 4.5. for the time 30.0. at each 
Reynolds number. There are two figures shown in each Fig. in which both figures 
are similar but second one is the in large version of first. 

The main problem in this method we faced that, there are lot of fluctuations 
in the calculation of pressure distribution due to the numerical noise. We tried a 
lot to get better results for pressure distribution but due to the numerical noise we 
could not get. The fairly large time step was a reduction in numerical noise. But by 
this idea, we have developed a general procedure by which we can obtain the flow 
diagram of any two dimensional multiple bodies, in which the CPU time is very less. 



67 ' 


Tmax=30 



Fig. 4.2. Full vortex cloud solution for two multiple bodies at Re - 20000 
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Fig. 4.3. Full vortex cloud solution for two multiple bodies at Re 


30000 
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Tmax=30 
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Fig. 4.4. Full vortex cloud solution for two multiple bodies at Re = 40000 
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Tmax=30 




Fig. 4.5. Full vortex cloud solution for two multiple bodies at Re = 50000 



Chapter 5 


CONCLUSIONS 

In the first part of the present work, we developed a procedure for the simulation 
of viscous diffusion in a boundary-layer flow. An operator splitting algorithm has 
been developed in which the viscous diffusion and convection are treated one after 
the other. The random walk technique is used for simulation of viscous diffusion. 
The boundary-layer velocity profile is calculated at one station and compared with 
theoretically obtained profile. At these values of moderate Reynolds numbers (500, 
1000 and 2000), the computed velocity profile is in fairly good agreement with 
the Blasius profile. The velocity profile has also been calculated at large Reynolds 
numbers (50000, 100000), for boundary-layer on a flat plate. It was observed that 
the velocity profile obtained is typical of the turbulent boundary-layer, and there 
appears to be a fair agreement with the law of the wall. 

The second part of the thesis is the simulation of the flow past a bluff body 
using the full-vortex cloud modelling by the surface vorticity boundary integral 
method. It develops a discrete vortex element model wherein the body is represented 
by a number of discrete vortex panels, and vortices in the wake are discretized 
into individual vortices being shed periodically. The strengths of the vortices are 
related directly to the tangential velocity outside the boundary-layer.- The tangential 
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boundary condition is satisfied at a number of collocation points through an integral 
equation. The full-vortex cloud model handles both diffusion and convection of 
vortices from the surface of the body. The convection of all discrete vortices takes 
place under the influence of all other discrete vortices and uniform stream velocity, 
whereas the viscous diffusion of all shed discrete vortices is simulated by using 
random walk technique. In this model vortices are continuously shed from the 
surface of the body into the flow fields and are free to move under the combined 
influence of vortex cloud and the external flow field. Finally the drag and lift 
coefficients are calculated by getting pressure distribution on the surface of the body 
and compared with the available data. We get good representation of the nature of 
the wake but the pressure calculation suffer from too much numerical noise. 

In the last part of this thesis an attempt is made to develop a procedure for 
calculating the flow past multiple two-dimensional bodies. We assembled two mutu- 
ally interacting bodies and calculated the combined effect after carefully considering 
the mutual interaction. We got the satisfactory vortex cloud solution for various 
Reynolds numbers. The pressure distribution and force coefficients are not satisfac- 
tory 

It appears that this formulation has the following advantages: 

(1) The CPU times involve are an order of magnitude long if we adopt any other 
method of simulation of viscous flow. 

(2) Change in shapes and conditions are relatively easily handled in this proce- 
dure so that quick calculation are possible with changed geometries. 

The main disadvantage in the method lies in the fact that calculation of pressure 
is still not fully understood. There is a lot of numerical noise in pressure calculation 
in unsteady flow, and the results are not yet satisfactory. 
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